Acoustic wave propagation in a macroscopically 
inhomogeneous porous medium saturated by a fluid 

L. De Ryck* J.-P. Grobyj P. Leclairef W. Lauriks§ 
A. Wirgin^C. Depollierll and Z. E. A. Fellah** 

February 2, 2008 



Abstract 

The equations of motion in a macroscopically inhomogeneous porous 
medium saturated by a fluid are derived. As a first verification of the 
validity of these equations, a two-layer rigid frame porous system con- 
sidered as one single porous layer with a sudden change in physical 
properties is studied. The wave equation is derived and solved for this 
system. The reflection and transmission coefficients are calculated nu- 
merically using a wave splitting-Green's function approach (WS-GF). 
The reflected and transmitted wave time histories are also simulated. 
Experimental results obtained for materials saturated by air are com- 
pared to the results given by this approach and to those of the classical 
transfer matrix method (TMM). 

Introduction. — An inhomogeneous medium is one with properties that vary 
with position. Inhomogeneous and layered materials can be found in many 
fields of physics, for instance in optics and electromagnetism El E] or in 
acoustics 13 12] . Other examples are geophysical media jH] , granular [Tj [S] or 
porous materials with depth-varying physical properties. The study of the 
acoustic wave propagation in inhomogeneous porous and granular media is 
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of great interest in building engineering, in petroleum prospection and civil 
engineering. 

Acoustic wave propagation in fluid-saturated porous media is relatively 
well-known thanks to the early work of Biot ^ and to the contribution of 
many researchers [lOl 111! I12j since then. A porous medium can be defined 
as a biphasic material with a solid phase and a fluid phase. This definition 
encompasses sandstones, concrete, absorbing polyurethane foams, snow or 
bones for example. While homogenized porous media have been extensively 
studied, macroscopically-inhomogeneous porous media have received far less 
attention. 

In this article, acoustic propagation in fluid-saturated macroscopically- 
inhomogeneous porous materials is studied. It is assumed that the wave- 
lengths are greater than the average heterogeneity size at the pore scale 
so that the physical properties of the material are homogenized. However, 
the homogenized properties can vary with the observation point within the 
material at the macroscopic scale of the specimen. The equations of mo- 
tion are derived from Biot's alternative formulation of 1962 |13j in which 
the total stress tensor, the fluid pressure, the solid displacement and the 
fluid/solid relative displacement w are used. It was briefly stated by Biot 
jl3j and confirmed [T^ that these variables should be employed to describe 
the acoustical properties of porous media with inhomogeneous properties. 
Among many possible applications, this work is a first contribution towards 
the determination of the inhomogeneity profile of unknown materials by the 
use of inversion methods. 

The equations of motion for a macroscopically inhomogeneous medium 
with elastic skeletton are derived. A first verification of the validity of the 
proposed equations is the study of a porous material saturated by air in the 
rigid frame approximation ^21; i-e- when the fluid is light and the solid skele- 
ton therefore relatively immobile. The porous material is then considered 
as an equivalent fluid with frequency-dependent and depth-dependent effec- 
tive density and bulk modulus. In this case, a wave equation is derived and 
solved numerically for a two-layer porous system treated as one single porous 
medium with a sudden but continuous change in physical properties. This 
provides an excellent means of comparing the proposed method (the WS- 
GF method) to the results of the well established Transfer Matrix Method 
(TMM) developed to calculate the acoustical properties of multilayer porous 
systems The WS-GF method is applicable to any depth-dependent in- 
homogeneous system and the two-layer system is just chosen for a testing 
purpose. An inhomogeneity function with a shape similar to the Heavi- 
side function, but which is analytical, is used. The jump with controlable 
steepness is created by multiplying the parameters by the inhomogeneity 
function. The thickness of the inhomogeneous medium is equal to the sum 
of thicknesses of each layer. 

The reflected and transmitted pressure fields are calculated from a known 
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incident plane wave impinging at normal incidence using a Wave Splitting- 
Green's functions approach (WS-GF) 3 . The method can be extended 
to include oblique incidence by working on the appropriate component of 
the wavevectors. The reflection and transmission coefficients and the time 
histories of the reflected and transmitted waves simulated with this technique 
are compared to the results given by the classical transfer matrix method 
(TMM) for multilayer porous materials Experimental results obtained 
in a two-layer material saturated by air are compared to the simulations. 

The equations of motion in an inhomogeneous poroelastic medium. — The 
constitutive linear stress-strain relations in an initially stress- free, statistically- 
isotropic porous medium can be written as jl3j 

aij = 2fieij + 6ij{XcB-aMC), (1) 
P = Mi-a0 + C), (2) 

where aij is the total stress tensor and p the fluid pressure in the pores; 
Sij denotes the Kronecker symbol (the summation on repeated indices is 
implied); 9 = V ■ u and C = —V • w are respectively the dilatation of the 
solid and the variation of fluid content where u is the solid displacement 
and w = (U — u) the fluid/solid relative displacement (U is the fluid 
displacement); (j) is the porosity; Cij = \ {uij + uj^i) the strain tensor of the 
solid (the comma denotes spatial partial derivatives); Ac = A + a^M, where 
A, n, M are elastic constants and a a coefficient of elastic coupling. These 
parameters were defined by Biot and Willis jl5| . 

Applying the momentum conservation law in the absence of body forces, 
the equations of motion are written 

V-<T = pii + pfW, (3) 
—Vp = p/ii-l-mwH — Fw, (4) 

where the dot and double dots notations refer to first and second order time 
derivatives, respectively; /)/ is the density of the fluid in the (interconnected) 
pores, p the bulk density of the porous medium, such that p = {l — (f))ps+(f)pf 
where ps is the density of the solid; m = PfT^/(f) is a mass parameter 
defined by Biot ^Sl; ''"oo is the tortuosity, rj the viscosity of the fluid, n the 
permeability and F the viscosity correction function. This function has been 
studied in detail by Johnson et al. ^U] and by Allard jl2j . 

We shall now derive the equations of motion for an inhomogeneous 
porous layer or a half space whose properties vary along the depth x. The 
fact that the medium is inhomogeneous is not incompatible with the fact 
that it is isotropic as the inhomogeneous medium can be considered as a 
superposition of an infinite number of thin isotropic sub-layers of thickness 
dx. Therefore, the following parameters in the above equations are now 
dependent on x: A, Ac, a, M, 0, p, m. Too, k and F. The ratio rj/n is 
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the flow resistivity and is often used instead of k in engineering acoustics 
applications. It is denoted by Rf here. The viscosity correction function F 
incorporates the viscous characteristic length A of Johnson et al. [lOj and 
the thermal characteristic length A' of Champoux and Allard jl6j . These 
parameters also depend on x. Inserting equations Q and ((2j into equations 
(jSJ and Q) yields the equations of motion in term of the displacements 



' V[(Ac+2^)V-u + aMV-w]-VA[^VAu]- 
2V/uV-u+2V/iA(VAu) + 2 [V^-V] u = pu+p/w, 

V [MV-w + aMV-u] = p/ii + mw + -Fw, 



(5) 



where the x-dependence of the constitutive parameters has been removed to 
simplify the notations. 

Wave equation in a rigid frame porous medium. — The previous equations 
can be applied to porous media with an elastic frame. Under the assump- 
tion of a rigid frame, u = and equations can be simplified. The 
porous medium can be considered as an equivalent fluid (at the scale of the 
wavelengths) described in the frequency domain by 



jujp = Ke{x,u))V .[4>{x)\J], 

-Vp = jujpeix,Uj)(f){x)\J, 



(6) 
(7) 



where pe{x,u;) and K(,{x,uj) are respectively the effective density and bulk 
modulus of the inhomogeneous equivalent fluid. Their expressions are 



Pe{x,u:) = pf 



4>{x) 



1 



■J^^^Fix,.) 

tOPfToo{xj 



Ke{x,Uj) 



7^0 



7-(7-l) 



^G{x,B^uj] 



(8) 
(9) 



where 7 is the specific heat ratio, Pq the atmospheric pressure and the 
Prandtl number. F{x,uj) and G{x, B'^uj) are the well-defined correction 
functions of the Johnson- Allard model |1U1 112j . The effective density and 
bulk modulus of the inhomogeneous equivalent fluid are functions of the 
frequency-independent parameters (/'(x), roo(a;), Rf{x) = q/n^x), A{x) and 
A'(x). These are the parameters that should be multiplied by the inhomo- 
geneity function in order to account for the change in properties. 

The wave equation in p is obtained by combining equations and © 



uj'^p + Ke{x, u>)V 



Pe{x,Uj) 



Vp 



0. 



(10) 



Wave splitting technique. — The second order differential operator of the 
wave equation in an homogeneous fluid can be factorized and this yields a 
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system of two coupled first order differential equations 



dx T 



Co 



p = ±Const X — p 
Co 



(11) 



where cq is the sound speed in air, p"^ corresponds to right-going waves and 
p~ corresponds to left-going waves, the sum of which equals p. This is the 
so-cahed Wave Splitting description, which was mainly used in scattering 
problems in the time domain in electromagnetism and then adapted 
to the frequency domain [3]. It can be seen as a change of variables from 
{p,dxp) to 




p"(L,o)) = p^ 



Air-coupled 

ullrasonic 

transducer 



Figure 1: Slab of inhomogeneous porous material. 

An inhomogeneous porous slab on which impinges an incident wave is 
shown in Fig. ^ Applied to the wave equation (fT(l|) . the wave splitting 
transformation is 

' ' p ± Zoct){x)\J .n (12) 



1 r 



where Zq = pjCQ is the characteristic impedance of the fluid surrounding the 
slab and n the unit normal vector (Fig. A system of linear first order 
coupled differential equation is obtained by combining equation ()12() with 
equations and ©: 

i dxP^ = A+{x,iv)p+ + A-{x,uj)p-, 
I dxP~ = —A~{x,u!)p'^ — A'^{x,u>)p~ 
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with 



zb 



(14) 



Numerical resolution of the wave equation. — The computation principle 
is the following: are first calculated in the surrounding homogeneous 
fluid at X = L. An infinitely thin homogeneous layer of thickness dx is 
then inserted at x = L — dx. The characteristic impedance of this layer 
is the same as the characteristic impedance of the inhomogeneous material 
at that point. At x = L a new set of are determined with the help of 
equations H13() . A new thin homogeneous layer is added at x = L — 2dx with 
the corresponding values of pe and K^. Using the updated values of p^ at 
X = L, the pressure subfields p^{L — dx,uj) are calculated. The operation 
is repeated until the last infinitely thin layer is added at a; = 0. For each 
addition of a new layer, the continuity conditions on p and (/)(j;)U • n are 
implicitely accounted for on both sides of the cumulated slab. 

Green's function approach. — The initialization of the procedure described 
above requires that p^ must be determined at x = L. To avoid this calcula- 
tion, a Green's function approach |17|in] can be used. Two Green's functions 

are defined by 

Vxe[0,L], p^{x,uj) = G'^{x,u;)p'^{L,uj). (15) 

Green's functions are characteristic of the sole material properties and de- 
scribe the internal field within the material. The boundary conditions at 
X = L are known and are G~^{L,ij) = 1 and G^{L,uj) = 0. The system of 
coupled first order linear differential equations in G^ obtained by inserting 
in l^l'A^ can be solved numerically using a Runge-Kutta routine. 
The reflection and transmission coefficients R{lij) and T{uj) are deduced 
from p^ 

p-{0,u;) = R{lo)p+{0,lo), (16) 
p+(L,c^) = T{Lo)p+{0,io). (17) 

From these coefficients, the reflected and transmitted waves can be simu- 
lated. 

In the numerical simulations, the function chosen to create the changes in 
physical properties is a distribution obtained from integrating the Gaussian 
normal distribution: I{x)=C{l — erf{—{x — xo)/r)) where C is a constant, 
xq the position of the jump and r a steepness factor. The steeper is the jump 
the finer the stepping dx must be for better accuracy. In the simulations, 
400 points were chosen to discretize the total slab and dx = 17.1/400 = 
0.0428mm. The value chosen for r was r = O.ldx. Values of r less than 
this value had little effect on the computed results. Smoothing the jump by 
taking r = Wdx resulted in an important reduction of the signal reflected 
at the interface between the two layers. 
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Experimental results - comparison with predictions. — The experimental 
principle is shown in Fig. ^ where an airborne ultrasonic wave is gen- 
erated and detected at normal or oblique incidence by specially designed 
(ULTRAN) air-coupled transducers in a frequency range between 150 and 
250kHz. The incident wave is partly reflected, partly transmitted and partly 
absorbed in the inhomogeneous porous layer. The materials studied are 
highly porous polyurethane foams saturated by air. The layers are put in 
contact, not glued. The physical parameters used for each layer were deter- 
mined using a previous method 18 and are displayed in Table I. 










A 


A' 




Thickness 








(^m) 


(fim) 


(Ns.m-^) 


(mm) 


Layer 1 


0.96 


1.07 


273 


672 


2843 


7.1 


Layer 2 


0.99 


1.001 


230 


250 


12000 


10.0 



Table 1: Properties of the two-layer medium studied. 



The reflection experiment was carried out with a single transducer used in 
the pulse-echo mode |19j . Another transducer was required on the other side 
of the specimen to measure the transmission. A particularly good agreement 
between the experimental results and the results of the simulations is found 
for the reflected and transmitted signals (Fig. I2ta) and IHb)). In Fig. 
|2fb), the three curves are almost indistinguishable. The waveforms were 
calculated from the experimental incident waveform, which was recorded and 
introduced into the simulation routines. The agreement is also very good for 
the reflection and transmission coefficients (Fig. ^ in the frequency range 
170 — 230kHz. In all figures, the WS-GF-curves cannot be distinguished 
from the TMM-curves. The discrepancies below 170kHz and above 230kHz 
in Fig. 131 are attributed to the limited bandwidth of the incident signal, 
making the experimental results more sensitive to noise outside the useful 
bandwidth. 
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Figure 2: Experimental signal (a) reflected and (h) transmitted by a two- 
layer porous system. Comparison with the signals simulated by the WS- GF 
method and by the TMM method. The incidence is normal. 
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Figure 3: Modulus of the experimental normal incidence reflection and trans- 
mission coefficients R{uj) and T{uj). Comparison with the reflection and 
transmission coefficients simulated by the WS-GF method and by the TMM 
method. 

Conclusion. — The equations of motion in fluid-saturated inhomogeneous 
porous media were derived from Blot's second formulation of 1962 |13j . A 
wave equation, valid in the rigid frame approximation, was also proposed 
and solved for a two-layer system of homogeneous porous material consid- 
ered as one single inhomogeneous porous layer with a sudden change in its 
physical properties. The transition from the properties of the first layer to 
those of the second layer was modeled by continuous inhomogeneity function 
with the shape of the Heaviside function. Excellent agreement was found 
between experimental and simulated waveforms and also between their re- 
spective reflection and transmission coefficients in a frequency range between 
170 and 230kHz. Future developments of this work are oblique incidence 
experiments and the study of other types of inhomogeneity profiles such as 
the linear profile. This work should contribute to the determination of the 
inhomogeneity profile of unknown materials by the use of inversion methods. 
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